Pair diffusion, hydrodynamic interactions, and available volume in dense fluids 
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We calculate the pair diffusion coefficient D{r) as a function of the distance r between two 
hard-sphere particles in a dense monodisperse suspension. The distance-dependent pair diffusion 
coefflcient describes the hydrodynamic interactions between particles in a fluid that are central to 
theories of polymer and colloid dynamics. We determine D{r) from the propagators (Green's func- 
tions) of particle pairs obtained from discontinuous molecular dynamics simulations. At distances 
exceeding ~3 molecular diameters, the calculated pair diffusion coefficients are in excellent agree- 
ment with predictions from exact macroscopic hydrodynamic theory for large Brownian particles 
suspended in a solvent bath, as well as the Oseen approximation. However, the asymptotic 1 /r dis- 
tance dependence of D{r) associated with hydrodynamic effects emerges only after the pair distance 
dynamics has been followed for relatively long times, indicating non-negligible memory effects in the 
pair diffusion at short times. Deviations of the calculated D{r) from the hydrodynamic models at 
short distances r reflect the underlying many-body fluid structure, and are found to be correlated 
to differences in the local available volume. The procedure used here to determine the pair diffusion 
coefflcients can also be used for single-particle diffusion in conflnement with spherical symmetry. 



I. INTRODUCTION 

Pair diffusion features prominently in theories of 
reaction-diffusion dynamicfP^ describing processes where 
reactant encounters are required, such as ligand binding 
and aggregation or fluorescence quenching. The hydrody- 
namic interactions quantified by the distance-dependent 
diffusion coefficient are also central to the theory and sim- 
ulation of polymer dynamics, including protein folding 
simulations in implicit solvent, the hydrodynamic cou- 
pling in dense colloidal suspensions, and the function 
of nanomachines and bacterial flagella.^ Considering the 
broad importance of pair diffusion in theories of molec- 
ular kinetics, it may seem surprising that little is known 
about the pair diffusion coefficient and its dependence 
on the particle distance. Formidable challenges in both 
theory and simulation^^HSl have resulted in often contra- 
dictory results for this fundamental quantity. 

Theoretically, the pair diffusion coefficient D{r) (with 
r the distance between two particles) has been at- 
tacked from two opposite directions, building up from 
kinetic theory'^' or projecting down from macroscopic 
hydrodynamics .'^CJ For D{r), kinetic theory had limited 
success at high fluid packing densities, largely because of 
the complexity of the molecular motions in dense fluids 
resulting from their many-body character. At the other 
extreme, details of the molecular structure of the sol- 
vent are ignored in estimates of the pair friction derived 
from macroscopic hydrodynamics, for instance by using 
the Oseen or Rotne-Prager tensors.^^El Nevertheless, this 
approach has proved useful in studies of the dynamics 
of large and sufficiently distant pairs of colloidal parti- 
cles in a solvent,'^ where macroscopic hydrodynamics is 
expected to apply; but it is not immediately applicable 
when solute and solvent particles are of comparable size, 
for instance in (aqueous) solutions of (bio)polymers. 

Here, we determine the pair diffusion coefficient di- 



rectly from the simulated many-body dynamics in a 
dense fluid. We focus on particles of the same size as 
the solvent molecules. This small-solute regime is of par- 
ticular relevance because, on the one hand, it allows us to 
quantify hydrodynamic interactions relevant for molecu- 
lar motions, including the dynamics of (bio)polymers in 
solution, and, on the other hand, it is far outside the 
regime where macroscopic hydrodynamics should be ex- 
pected to apply. 

The paper is outlined as follows. In section I, we de- 
scribe the methodological details, including the theory 
to calculate the pair diffusion tensor, the algorithm used 
to determine the required Green's functions from simula- 
tion data, the simulation parameters, and the validation 
procedure. We validate our method by computing the 
pair diffusion coefficient for two spherical particles sub- 
ject to Brownian dynamics. In the results section II, 
we first present a comparison of Green's functions ob- 
tained from simulations against those predicted from our 
diffusion model, finding excellent agreement over 8 or- 
ders of magnitude. Then we examine the pair diffusion 
coefficients as a function of distance between two parti- 
cles for several fluid packing fractions, and compare the 
simulation results to the predictions of macroscopic hy- 
drodynamic theories. Finally, we show that the position- 
dependent pair diffusion coefficient is correlated to the 
local available volume. In the Appendix, we discuss the 
calculation of the angular pair diffusion coefficient. 



II. METHODS 
A. Theory 

In the following we present the theory to calculate the 
position-dependent pair diffusion tensor from simulation 
trajectory data. The diffusion tensor D of the vector r 
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between two given particles in an isotropic and homoge- 
neous fluid has spherical symmetry: 



D(r) = i:'_L(r-)e^e^ + D« (r) (e^ee 



(1) 



where r = |r| is the length of the pair vector; D±{r) and 
D||(r) are the scalar diffusion coefficients in the radial 
and tangential directions, respectively; and e^, e^, and 
are the orthonormal unit vectors of the spherical polar 
coordinate system, with pointing in the radial direc- 
tion, and eg and being tangential to longitudes and 
latitudes, respectively. The Smoluchowski (or Fokker- 
Planck) equation describing the diffusion of the pair vec- 
tor then takes on the following form: 



dtp - div [D(r)e-'^^grad {e^^ p)] 



(2) 



where p = p{r,9,t\r' ,6q = 0,i = 0) is the Green's func- 
tion for a pair vector starting at a distance r' and az- 
imuthal angle 6*0 = 0, without loss of generality because 
of the isotropic space (making the ip distribution uni- 
form); V{r) is the distance-dependent free energy surface; 
P — (kBT)^^ is the inverse temperature; dt is the partial 
derivative with respect to time; and "div" and "grad" are 
the divergence and gradient operators in spherical polar 
coordinates, respectively. 

We will in the following use x = cos 9 instead of 9. Let 
P{r, X, t\r' , 0) be the Green's function in terms of this new 
variable. The diffusion equation Eq. ^ then becomes: 



dtP^dr[D^{r) {13V' + dr 

+^d4{i-x')d^p] 



p] 



(3) 



where V — dV{r)/dr. By integrating over x — cos 9 we 
obtain a diffusion equation for the Green's function in 
the radial direction alone, with the second term on the 
right hand side vanishing: 



dtG = dr [D±{r) {l3V'G + drG)] 



(4) 



where G{r,t\r' ,0) — J^-^dx P{r,x,t\r' ,0) is the proba- 
bility for the pair distance to be in (r, r + dr) at time 
t, starting from r' at time 0. As a consequence, we 
can treat radial diffusion separately using standard one- 
dimensional diffusion, irrespective of the angular motion. 
In an appendix, we outline an extension of the theory to 
the orientational diffusion of the pair distance vector. 



Algorithm to determine pair distance diffusion 
coefficient 



Here we focus on the calculation of the position- 
dependence of the pair-distance diffusion coefficient 
D{r) = D±(r), where we have dropped the subscript 
for notational simplicity. In our calculations of D{r), we 
face the dual challenges that it depends on the particle 
distance r, and that the pair dynamics becomes diffusive 
only at times at which the influence is felt of the un- 
derlying free energy surface (or potential of mean force) , 
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FIG. 1. Pair diffusion coefficient D{r) for two freely diffus- 
ing Brownian particles of diameter 1 with periodic boundary 
conditions. Results for different grid sizes Ar = 0.004 (plus), 
0.032 (cross), 0.024 (star), 0.016 (square) are compared to 
the exact value 2Do — 0.1 (horizontal line). The vertical 
solid line marks the contact distance r = 1. To assess arti- 
facts from periodic boundary conditions, the vertical dashed 
lines mark distances r = L/2, L/\/2, and \/3I//2, where cen- 
tered spheres touch the faces, edges, and corners of the cubic 
simulation box, respectively. 



F{r) = -kB.T\TLg{r) = V[r) + 2kBT\nr, where g{r) is 
the pair correlation function of the two particles in the 
fluid. To disentangle the diffusive spread of the pair dis- 
tance distribution from the drift of the mean position 
as a result of the underlying free energy surface, we use 
the propagator (or Green's function) G(r, t\r\ 0)dr. In 
constructing a diffusion model, we assume that G satis- 
fies the Smoluchowski diffusion equation Eq. Q, where 
the term within the brackets is the negative of the radial 
probability flux. A spatial discretization of the Smolu- 
chowski equation^ results in a master equation that de- 
scribes the pair dynamics between neighboring intervals 
along r. The particle-pair trajectories in the simulations 
are discretized by assigning pair distances into bins i 
along r, and then counting the numbers Nji that a pair 
distance is in bin i at time r, and in bin j at time r + At, 
irrespective of its location at intervening times, with At 



the lag time. Nji is symmetrized, A^^^- 



Nji, consistent 



with microscopic time reversibility. We then find the pair 
diffusion coefficient D{r) that maximizes the path action 
of the observed discretized path. For the discretized dif- 
fusion model with given D(ri) and F{ri), the path ac- 
tion (or likelihood) L can be be written as a product of 
Green's functions that are expressed in terms of a ma- 
trix exponential.'^^ To optimize the action and find the 
diffusion model most consistent with the observed Nji, 
wc infer D{ri) and F{ri) using a Bayesian approachjIM 
with uniform priors in \nD{ri) and F{ri) ensuring scale 
invariance in time and space. 

In free diffusion, one typically fits a + QD^t (or, equiva- 
lently, 6I?o(^ + ''")) to the mean-square displacement, with 
the constant a (or the time shift r ~ a/GDo) accounting 
for initial fast molecular motions. Here, we employ a 
similar procedure by optimizing also the time origin r 
for transition counts Nij collected at several different lag 
times At, 2At, . . . , kAt = t, where t defines the "obser- 
vation time." 

To validate the procedure, we first run Brownian dy- 
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namics simulations for two spherical particles of unit di- 
ameter freely diffusing with diffusion coefficient Dq = 
0.05 in a cubic box of length L = 12.5 under peri- 
odic boundary conditions and with reflecting bound- 
aries at particle contact. By construction, in this case 
D{r) = 2Do, which is indeed recovered by the procedure 
for distances r < L/2 (Fig. [T]), nearly independent of grid 
size Ar. However, for r > L/2 and long lag times, the 
periodic boundary conditions cause artifacts because in 
the corners of the cubic simulation box the pair dynam- 
ics projected onto the minimum image distance depends 
not only on the length of the pair vector but also on its 
direction. 



III. SIMULATIONS 

To calculate D{r) for a particle pair in a dense fluid, we 
perform discontinuous molecular dynamics (DMD) sim- 
ulations of hard sphere (HS) fluids. In DMD, particles 
follow linear trajectories between collisions. In a colli- 
sion, the velocities of colliding particles are changed to 
conserve energy and momentum.^' To simplify the nota- 
tion, dimensionless quantities will be used, obtained by 
appropriate combinations of a characteristic length (HS 
particle diameter a) and time scale (a^JmP, where m 
is the particle mass). The packing fraction cf) = tt/j/G 
is defined in terms of the particle density p. To con- 
struct the Green's functions, we performed DMD sim- 
ulations with N = 2000 identical HS particles. Peri- 
odic boundary conditions were applied in all directions. 
The average self-diffusivity Dq was obtained by fitting 
the long-time {t ^ 1) behavior of the average mean- 
squared displacements Ar^ of the particles to the Ein- 
stein relation (Ar^) = 6Dot. To minimize the system- 
size dependence trajectories from simulations with 
N = 10000 particles were used to determine Dq, with 
remaining finite-size corrections of wl 

IV. RESULTS 

To test the applicability of the diffusion model, we 
compare its prediction for the dynamics of the pair dis- 
tance to actual simulation data collected over a range of 
time scales. Figure [2] shows that diffusion quantitatively 
captures the pair dynamics in the fluid. The Green's 
functions G{r,t\r' ,0) from the diffusion model and the 
results of the DMD simulation data are found to agree 
over 8 orders of magnitude. At the shortest observation 
time < = 1, we find that the Green's functions are es- 
sentially Gaussian with position-dependent widths. At 
longer times, t — 10 and 20, the underlying free energy 
surface shows its influence, distorting the propagators 
away from the Gaussian form expected for free diffusion 
on a flat surface. 

In Figure |3j we explore the effects of the spatial grid 
size Ar and the observation time t on the calculated pair 
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FIG. 2. Green's functions G{r,t\r' ,0) from simulations (sym- 
bols) and diffusion model (lines). G(r,t\r' ,0) is shown as a 
function of the pair distance r at packing fraction cj) = 0.325 
for time t = 1 (top panel), 10 (middle panel), 20 (bottom 
panel). We use an observation time of i = 20 to obtain 
diffusion model parameters, combining results for lag times 
At = 1, 2, ... , 20. The arrow in the top panel reflects increas- 
ing r' = 1,2,...,7. 



diffusion coefficient. We find that for Ar < 0.1, grid 
size effects are negligible. Figure [s] (bottom) shows that 
the effect of changing the observation time t is negligible 
only for shorter distances r < 3. In contrast, for longer 
distances D{r) is almost flat at a short observation time 
i = 4 and does not show the asymptotic l/r dependence 
expected from macroscopic hydrodynamic theory. How- 
ever, the expected 1 /r dependence is recovered for longer 
times t. This result implies that the hydrodynamic cou- 
pling at large distances is not instantaneous, such that a 
more accurate diffusion model would require the inclusion 
of memory effects in a frequency and position-dependent 
diffusion coefficient .^Sl For t > 16 the predictions are es- 
sentially independent of t. In all following calculations, 
we thus use Ar = 0.1 and t = 20. 

Having validated the procedure and diffusion model, 
we now examine the distance-dependent pair diffusion 
coefficients D{r) for different packing fractions (p. Fig- 
ure |4] (top panel) shows D{r) for the HS fiuid over 
a packing fraction range (f> = 0.325 — 0.48 (symbols 
from top to bottom). Also shown are the predictions 
for D{r) from hydrodynamic theory f or two spherical 
particles with slip boundary conditions'^ as well as 
the widely-used Oseen tensor correctiorP (for (f> — 0.4; 
dashed line), which for the pair diffusion coefficient is 
D{r) ~ 2Do — k-gT / {2TTr]r) where 77 is the solvent shear 
viscosity, taken from Ref. [13l We find that both the 
exact hydrodynamic theory and the Oseen approxima- 
tion (and similarly the Rotne-Prager tensor;- not shown) 
are remarkably accurate and quantitatively reproduce 
the large-r behavior. However, hydrodynamic predic- 
tions only qualitatively reproduce the observed decrease 
in D(r) near contact (r = 1) and lack any structure due 
to molecular correlations in the first- and second-shell 
around a particle. 

To characterize the effects of the molecular pack- 
ing structure on the pair dynamics, we plot in Fig. [4] 
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FIG. 3. Dependence of D{r) on diffusion model parameters. 
Pair diffusion coefficient D{r) versus distance r for a fiard- 
sphere fluid at packing fraction (j) = 0.35 obtained for different 
(top) grid sizes Ar (with fixed observation time t = 20) and 
(bottom) observation times t (with fixed grid size Ar — 0.1). 
The lag time is At = f consistently. 

(bottom panel) the normalized pair diffusion coefficient 
D{r)/2DQ for different packing fractions cj). As expected 
from macroscopic hydrodynamics, at large distances r 
the D{r)/2Do data collapse onto a single curve that is 
well represented by the hydrodynamic theory. Two im- 
portant observations are: (i) D{r)/2DQ is always less 
than 1, with pair diffusion slowed down by "hydrody- 
namic interactions." (ii) D(r) rises sharply just out- 
side distances of 1 and 2 particle diameters, and drops 
sharply just outside r = 1.5, and 2.5. This strong po- 
sition dependence, together with the short-time propa- 
gator G(r, i|r', 0) for Brownian dynamics being Gaussian 
with mean r = r' + t[D{r)f3drFr + drD], implies that at 
short times particle pairs just outside the first and sec- 
ond shell boundary "drift" outward, whereas those inside 
the boundary drift inward, beyond what is expected from 
the free energy gradient drF alone. The additional drift 
terms arise from the large gradients in D{r). We can 
understand this dynamic behavior (which does not vio- 
late microscopic time reversibility and detailed balance!) 
from the many-body packing effects. At r > 2, for in- 
stance, the interstitial space between the two particles 
is likely filled by a third one, which tends to drive the 
pair apart. In contrast, at r ^ 2, the interstitial space 
between the two particles is empty, and the two particles 
tend to move closer together. 

To gain further insight into the observed structure in 
D{r) and its relation to the static structure of the fluid, 
we plot in Figure [5^ both D{r) and the pair correlation 
function g{r). We find that there is some correlation 
between the structure in D{r) and g{r) except near the 
contact distance at r = 1 where these quantities are actu- 
ally anti-correlated. Somewhat counter-intuitively, this 
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FIG. 4. Pair diffusion for a hard-sphere fluid. (Top) Cal- 
culated pair diffusion coefficient D{r) versus distance r with 
increasing packing fraction — 0.325, 0.35, 0.375, 0.40, 0.42, 
0.44, 0.46, 0.48 (symbols, from top to bottom). Lines are 
the predictions of hydrodynamic theory (see text). (Bottom) 
Normalized pair diffusion coefficient D(r) /2Dq, where Dq is 
the self-diffusivity for a given 0. Symbols are ou r calcu lations, 
the thick line is the exact hydrodynamic theory,'^^'^ and the 
dashed line is the Oseen approximation. 



mostly positive correlation means that the pair diffusion 
is actually higher in the more densely packed regions. 
Similar behavior was observed for a HS fluid confined 
between hard walls where the local density was found to 
be strongly correlated with the local diffusion coefficient 
except near the walls."^ This behavior was found to be 
related to the physics of layer formation, with the avail- 
able volume, as probed by the local test-particle insertion 
probability Pqj being largest in the locally dense regions 
of space.!-' \ similar argument should hold in our case of 
a bulk HS fluid in which purely entropic excluded volume 
forces give rise to a structured g{r) profile to maximize 
the system entropy. The local insertion probability is 
given by PqIt) — p{r)/^ ~ pgir)/^, where the activity 
^ = exp(/3/i) / \^ is spatially invariant for an equilibrium 
fluid, with /i the chemical potential, and A the thermal 
wavelength. 

To test if D{r) is indeed related to P^ir), we calculate 
^ for the different packing fractions by utilizing grand 
canonical transition-matrix Monte Carlo simulations.!^ 
Figure [5}d shows D{r) versus Po{r) for different 4>. We 
find that the D{r) data approximately collapse onto a 
curve similar to the average bulk relationship (2_Do versus 
Po) that ignores any r dependence. Therefore, at least 
as a rough approximation, the local available volume can 
describe the pair diffusion in this case. 
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FIG. 5. Relation between fluid structure and dynam- 

ics. (Top) Pair diffusion coefficient -D(r) (symbols connected 
by lines) and scaled pair correlation function g' (r) = g{r-)/a 
(lines) versus distance r where a is an arbitrary scaling factor 
used to match D{r) and g' {r) at large r {(p = 0.325, 0.375, 
0.42, 0.48 from top to bottom). (Bottom) D{r) as a func- 
tion of the local fractional available volume Po{r) (symbols; 
increasing packing fractions from right to left). The line is 
2Do versus Po averaged over the entire system. 



V. CONCLUDING REMARKS 

The results of this paper shed Ught on the microscopic 
origins of the distance dependence of hydrodynamic in- 
teractions, in particular the role of particle packing and 
many-body motions, and help establish a range of vaUd- 
ity for the assumption of macroscopic hydrodynamics in 
the modeling of processes ranging from polymer dynam- 
ics to nanomachines, colloidal dynamics, and bacterial 
swimming. In practical applications, such as the calcu- 
lation of diffusional encounter rates, the significant devi- 
ations between the calculated pair diffusion coefficients 
D{r) and the ideal (and widely used!) assumption of 
D{r) = 2Dq = const, can result in substantial errors, 
with D{r) < 2Do consistently. At the least one should 
use a hydrodynamic theory, with both the exact theory 
and the Oseen tensor giving remarkably accurate results 
for hydrodynamic interactions at larger distances, and 
rough approximations in the regime dominated by molec- 
ular packing near contact. 



Appendix: Angular diffusion coefficient 

To treat the angular diffusion of pair distance vectors 
(or other vectors in an isotropic space), we notice that 
the second term on the right hand side of Eq. ([s]) corre- 
sponds to the angular momentum operator in quantum 
mechanics. We thus make the ansatz P{r, x,t\r' ,0) = 



S/^o C^iPi{^)qi{^: 0)i where the Pi{x) are the Legen- 
dre polynomials of order I, and the coefficients Ci do not 
depend on t and r. With this ansatz, we obtain uncou- 
pled one-dimensional evolution equations for each of the 
qi (with; = 0,1,...): 



dtqi 



= dr [D^{r) i/3V'qi + drqi)]-^^4^lil+l)qi . (A.l) 

For I — 0, this expression is identical to Eq. Q; for Z > 0, 
these are sink (or birth-death) equations for the qi , with 
sink terms whose strength increases quadratically with 
I, and with £'||(r)/r^. That is, at long times only the 
distribution uniform in x survives (with Po{x) — 1). 

Expressed in terms of Dirac (5-functions, the initial con- 
dition for the Green's function is P{r,x,t = 0|r'',0) = 
6{r — r')S{l — x) (where we chose the coordinate sys- 
tem such that the polar axis points in the direction 
of the pair distance vector at time zero), with nor- 
malization J_^dx J dr P(r, X, i|r', 0) = 1. By using 
the orthogonahty relations of the Legendre polynomi- 
als, J^^dxPi{x)Pmix) = 26im/i'2.l + 1) with 5im the 
Kronecker-5, we obtain 



P(x,r,t|r',0) 



21 + 1 



1=0 



Pi{x)qi{r,t\r',0) (A.2) 



where the qi satisfy Eq. (A.l I with initial conditions 
qi{r,0\r',0)=S{r-r'). 

For the sake of completeness, we also sketch an algo- 
rithm to obtain the distance-dependent radial and an- 
gular diffusion coefficients D±{r) and D^\{r) from sim- 
ulation data (or, equivalcntly, from experimental data, 
such as those obtained in colloidal-particle tracking ex- 
periments). 

1. Use counts of transitions Nj^ from bins i to j in 
the radial direction only (irrespective of the angu- 
lar motion) as input in the algorithrrP^l described 
above to calculate the one-dimensional position- 
dependent diffusion coefficients D± (r) , and the po- 
tential of mean force V{r). 

2. Determine counts Nja^i for transitions from bin i 
in the radial direction to bin j,a in a two dimen- 
sional histogram. Radial bins are indexed by j, 
and angular bins by a according to the cosine of 
the azimuthal angle, 



x{t) = cos 61 (t) 



r{t) ■ r(0) 

m\m\ 



(A.3) 



(with 6(0) = and x{0) = 1 by definition of the 
coordinate system). 

3. With D±{r) and V{r) already determined in the 
first step, the Green's function Eq. (A.2 1 can be 
calculated for a given estimate of -Dy (r) from a spa- 
tially discretized version^ of the sink equations, Eq. 
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(A.ll. With this Green's function, one can again 



use a Bayesian inference procedure (or maximum- 
hkelihood method) to estimate the D^^{r) (on lat- 
tice points halfway between the bin centers) that is 
most consistent with the observed transition counts 



Note that the infinite sum over I in Eq. (A. 2 1 has to 



be truncated in practical calculations. Note further that 
the same algorithm can also be used to determine the 
diffusion coefficients of a single particle in confinement 
with spherical symmetry. 
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